function sift_collect = SIFT
    sift_collect.init_img        = @init_img
    sift_collect.build_gauss_pyr = @build_gauss_pyr;
    sift_collect.build_dog_pyr   = @build_dog_pyr;
    sift_collect.scale_space_extrema = @scale_space_extrema
end

% create init img
function dbl_img = init_img(pic, img_dbl, sigma)
    pic = rgb2gray(pic);
    if(img_dbl)
        sig_diff = sqrt(sigma * sigma - 0.5 * 0.5 * 4);
        w = round(sig_diff * 3);
        if mod(w,2) == 0
            w = w + 1;
        end
        dbl_img = imfilter(imresize(pic, 2), fspecial('gaussian', [w,w], sig_diff), 'replicate');
    else
        sig_diff = sqrt(sigma * sigma - 0.5 * 0.5);
        w = round(sig_diff * 3);
        if mod(w,2) == 0
            w = w + 1;
        end
        dbl_img = imfilter(pic, fspecial('gaussian', [w,w], sig_diff), 'replicate');        
    end
end

function gauss_pyr = build_gauss_pyr(pic, scls, octvs, sigma)
    % scls  : number of scale spaces
    % octvs : number of octaves per scale
    sz = size(pic);
    gauss_pyr = cell(scls, octvs + 3);
    
    % sigmas
    k = 2^(1 / octvs);
    sig = zeros(octvs + 3, 1);
    sig(1) = sigma;
    sig(2) = sigma * sqrt(k * k - 1);
    for i = 2:octvs+3
        sig(i) = sig(i-1) * k;
    end
    % build 
    for s = 1:scls
        for i = 1:octvs+3
            if s == 1 && i == 1
                gauss_pyr{s,i} = pic;
            elseif i == 1
                gauss_pyr{s,1} = imresize(gauss_pyr{s-1, octvs}, 0.5);
            else
                w = round(sig(i) * 3);
                if mod(w,2) == 0
                    w = w + 1;
                end
                gauss_pyr{s,i} = imfilter(gauss_pyr{s,i-1}, fspecial('gaussian', [w,w], sig(i)), 'replicate');
            end
        end
    end
end

function dog_pyr = build_dog_pyr(gauss_pyr)
    sz = size(gauss_pyr);
    dog_pyr = cell(sz(1),sz(2)-1);
    for s = 1:sz(1)
        for o = 1:sz(2)-1
            dog_pyr{s, o} = gauss_pyr{s, o+1} - gauss_pyr{s, o};
        end
    end
end



function ex_seq = scale_space_extrema(dog_pyr, scls, octvs, contr_thr, curv_thr)
    prelim_contr_thr = 0.5 * contr_thr / octvs * 255;
    for s = 1:scls
        for o = 2:octvs+1
            sz = size(dog_pyr{s,1});
            for r = 5:sz(1)-5
                for c = 5:sz(2)-5
                    % contrast threshold
                    if abs(dog_pyr{s,o}(r,c)) > prelim_contr_thr
                        if is_extremum(dog_pyr, s, o, r, c)
                            feat = interp_extrema(dog_pyr, s, o, r, c, octvs, contr_thr);
                            %is_bad_curv()
                        end
                    end
                end
            end
        end
    end
    ex_seq = [];
end

function flag = is_extremum(dog_pyr, scl, octv, r, c )
    val = dog_pyr{scl,octv}(r, c);
    % is maximum
    if val > 0 
        for i = -1:1
            for j = -1:1
                for k = -1:1
                    if val < dog_pyr{scl,octv+i}(r + j, c + k)
                        flag = 0;
                        return; 
                    end
                end
            end
        end
    % is minimum
    else
        for i = -1:1
            for j = -1:1
                for k = -1:1
                    if val > dog_pyr{scl,octv+i}(r + j, c + k )
                        flag = 0;
                        return; 
                    end
                end
            end
        end
    end
    flag = 1;
end

function feat = interp_extrema(dog_pyr, scl, octv, r, c, octvs, contr_thr)
    while i < 5
        x = interp_step(dog_pyr, scl, octv, r, c);
        if abs(x(1)) < 0.5 && abs(x(2)) < 0.5 && abs(x(3)) < 0.5
            break;
        end
        r = r + round(x(2));
        c = c + round(x(3));
        octv = octv + round(x(1));
        if octv < 1 || octv > octvs || c < 5 || r < 5  ||
            c >= size(dog_pyr{scl,0})(2) - 5  ||
            r >= size(dog_pyr{scl,0})(1) - 5
                feat = 0;
            return;
        end
        i = i+1;
    end
    if i >= 5
        feat = 0;
        return;
    end

    contr = interp_contr(dog_pyr, scl, octv, r, c, x);
    if abs(contr) < contr_thr / intvls
        feat = 0;
        return;
    end

    feat = zeros(2,1);
    feat(1) = (c + x(2)) * pow( 2.0, scl);
    feat(2) = (r + x(1)) * pow( 2.0, scl);
end

function interp_step(dog_pyr, scl, octv, r, c)

end